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ABSTRACT 

Context. Circumstellar disks are considered to be the birthplace of planets. Specific structures like spiral arms, gaps, and cavities are 
characteristic indicators of planet-disk interaction. Investigating these structures can provide insights into the growth of protoplanets 
and the physical properties of the disk. 

Aims. We investigate the feasibility of using molecular lines to trace planet-induced structures in circumstellar disks. 

Methods. Based on 3D hydrodynamic simulations of planet-disk interactions obtained with the PLUTO code, we perform self- 
consistent temperature calculations and produce N-LTE molecular line velocity-channel maps and spectra of these disks using our 
new N-LTE line radiative transfer code Mol3D. Subsequently, we simulate ALMA observations using the CASA simulator. We 
consider two nearly face-on inclinations, five disk masses, seven disk radii, and two different typical pre-main-sequence host stars 
(T Tauri, Herbig Ac) at a distance of 140 pc. We calculate up to 141 individual velocity-channel maps for five molecules/isotopoloques 
(12 c16o, HCO^, HCN, and CS) in a total of 32 rotational transitions to investigate the frequency dependence of the structures 

indicated above. 

Results. We find that the majority of protoplanetary disks in our parameter space could be detected in the molecular lines considered. 
However, unlike the continuum case, gap detection is not straightforward in lines. For example, gaps are not seen in symmetric 
rings but are masked by the pattern caused by the global (Keplerian) velocity field. By comparison with simulated observations of 
undisturbed disks we identify specific regions in the velocity-channel maps that are characteristic of planet-induced structures. 
Conclusions. Simulations of high angular resolution molecular line observations demonstrate the potential of ALMA to provide 
complementary information about the planet-disk interaction as compared to continuum observations. In particular, the detection of 
planet-induced gaps is possible under certain conditions. 
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1. Introduction and motivation 

It is commonly assumed that planets form in circumstellar disks 
around young pre-main-sequence (PMS) stars (e.g., |Mordasini| 
et al.|[20T0| . During this formation process planets perturb 
the disk becaus e of their gravitational potential ( [Goldreich &| 
Tremaine|1980| ). The results of this interaction are characteristic 
structures such as spiral arm s (|Ward||1997[ [Tanaka et al.||2002| ) 
and gaps ( jPapaloizou & Lin||1984| ). Therefore, observations of 
these large-scale structures provide indirect constraints on the 
existence of protoplanets and on the relevant physical processes 
involved in the protoplan et-d isk interaction. 

Wolf & D’ Angelo| ( |20Q5| ) and jRuge et al.| ( |2Q13] ) showed the pos¬ 
sibility of tracing gaps in continuum re-emission light for vari- 
ous disk configura tions and ALMA configurations. In addition, 
[Ruge et aT] ( |2014| ) investigated the feasibility of tracing gaps in 
continuum re-emission maps and scattered light images simulta¬ 
neously. In contrast to this, only a few studies exist investigat¬ 
ing the feasibility of tracing planet-induced structures through 
molecular line observations. [Semenov et al.| ( [2008] ) showed the 
feasibility of drawing conclusions about disk properties such 
as large-scale temperature gradients and c hemical stratification 
from simulated HCO^ (4-3) observations. Dutrey et"STl ( 2008 ) 
used spatially unresolved but spectroscopically resolved line ob¬ 


servations to trace inner cavities due to the radial dependence of 
the (Keplerian) velocity field. [Cleeves et al.[ ( [201 1] ) demonstrated 
the feasibility of detecting a huge inner cavity (R ~45 AU) us¬ 


ing molecular lines (CO, H2CO2). Using ALMA, [Casassus et al. 
( [2013| ) were able to image multiple molecular gas species within 
the dust free inner hole of HD 142527 and to utilize the vary¬ 
ing optical depths of the different species to detect a flow of gas 
crossing the gap opened by a planet. 

However, the main component of protoplanetary disks, cold 
molecular hydrogen (H2), is very difficult to observe directly 
in these disks because of the missing electrical dipole moment. 
Thus, other molecules with lower excitation temperatures must 
be observed to derive the dynamics, temperature, density, and 
chemical structure in these disks. The most frequently observed 
molecule is CO, mainly because of its high abundance and low 
excitation energy. Observations of various isotopologues of this 
molecule provided the first constraints on the physical properties 
such as the (vertical) gas temperature structure of young disks 
(e.g., [Pietu et al.|2007] ). Through comprehensive surveys of pro¬ 
toplanetary disks around solar-type stars other mol ecules such as 
HCO^, H2CO, CS, CN, and HCN were d etected jKastner et af 

2004| . Recently, Chapillon 


1997 


et al. 


Oberg et al. 


2010 


Thi et al. 


(2012a) provided evidence of HC3N in the disks around 
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GO Tau and MWC 480 using the IRAM 30 m telescope. How¬ 
ever, because of insufficient spatial resolution and the lack of 
sensitivity these observations were restricted to the simplest and 
most abundant molecules. The Atacama Large (Sub)-Millimeter 
Interferometer (ALMA) will potentially allow the direct study 
of the molecular composition and distribution in protoplane¬ 
tary disks in the submillimeter wavelength regime both at global 
scale as well as in the potential planet-forming region. The first 
ALMA observations of HD 163296 show that CO radiation orig¬ 
inates from a thin layer with an opening angle of about 15° with 
respect to the disk midplane ( [Rosenfeld et al.|20T3}|de Gregorio-| 
Monsalvo et al.|[MT3] ). More recently with the unprecedented 
sensitivity of ALMA, cyclic C 3 H 2 (which has been proposed as 
a useful probe of radiation penetration in disks) has been de¬ 
tected in the disk surrounding HD 163296 ( [Qi et al.|2013| ). 

The motivation to search for gaps in protoplanetary disks using 
molecular line observations is given not only by the aim to in¬ 
directly detect planets, but even more to derive complementary 
constraints on the disk physics from the planet-disk interaction. 
For example, the resulting gap structure, given by its shape and 
depth in the gas and dust distribution, depends not only on the 
mass of the planet and grain size, but also on the interaction of 
the planet with the gas and dust phase ( [Paardekooper & Mellema| 
|2004| ). As a result, larger dust grains (> 150 //m) are expected to 
decouple from the gas, pile up at the edges of the gap and signifi¬ 
cantly alter the gas-to-dust mass ratio ( [Paardekooper & Mellema| 
|2006| ). Moreover, we emphasize that a gap in the density dis¬ 
tribution of a cirumstellar disk is not necessarily an irrevocable 
sign of an embedded planet. [Flock et aT] ( |2015| ) showed that MRI 
turbulence is sufficient to generate a gap in the gas phase with- 
out the presence of a planet. Recently, [Bruderer et al.[ ( [2Q14[ ) 
were able to detect CO emission inside the dust cavity of IRS 
48. They concluded that a massive inner companion is the main 
driver for the clearing of the cavity and excluded photoevapora¬ 
tion and grain-growth only. Therefore, complementary observa¬ 
tions of the gas and the dust phase are crucial in order to improve 
and verify our understanding of the predominant physics; e.g., 
the interaction between the gas and dust component in the planet 
forming regions of protoplanetary disks. As a first step into this 
potentially new field of studying disk physics through planet- 
disk interaction, we investigate the feasibility of tracing planet- 
induced gaps through spatially and spectroscopically resolved 
molecular line observations. For this purpose we developed the 
N-LTE line and dust continuum radiative transfer code Mol3D, 
which is based on the radiative transfer code MC3D ([Wolf et al.[ 
[T999l[Wolfl2003l ). 


Our studies are based on spatial density and velocity distribu¬ 
tions resulting from 3D hydrodynamic simulations of planet-disk 
interaction using the PLUTO code ( [Mignone et al.[2007) ). The 
temperature structure and corresponding molecular line maps are 
calculated with MolSD. These velocity-chan nel maps are used as 
an input for the CASA (ver. 4.2) simulator ( [Petry & CASA De^ 
[velopment Team[[2012[ ). Finally, we investigate the feasibility 
of tracing planet-induced gaps in the simulated ALMA maps. 
Based on these results, our goal is to prepare future ALMA ob¬ 
servations of similar protoplanetary disks. In contrast to ear¬ 
lier studies, which were often focused on exemplary case studies 
(e.g., [Cleeves et al.][2Qll[ ), we explore a large parameter space 
covering various disk properties such as their total mass, size, 
and abundant molecules of typical pre-main-sequence T Tauri 
and Herbig Ae stars. 

This study is structured as follows. In Sect. we briefiy 
present the radiative transfer code MolSD, discussing its under¬ 
lying physical assumptions and test cases. In Sect, [^we describe 
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the study of tracing gaps using molecular lines with ALMA. 
Finally, in Sects. and we discuss our results and provide 
conclusions about the feasibility of tracing planet-induced gaps 
though molecular line observations with ALMA. 


2. Molecular line radiative transfer code 

In this section, we briefiy describe our line and dust continuum 
radiative transfer code MolSD. 


2.1. Mol3D features 


The N-LTE 3D parallelized line radiative transfer code Mol3D 
features full 3D spherical, cylindrical, and Cartesian grids with 
several types of boundary spacing (e.g., logarithmic, linear). 
MolSD typically uses three steps to produce velocity-channel 
maps, spectra, dust continuum maps, and spectral energy dis¬ 
tributions (SED). In the first step, the dust temperature is calcu¬ 
lated based on the assumption of local thermal equilibrium. To 
achieve a high efficiency, we use a combination of the continu¬ 
ous absorption method pr oposed by [Lucy \\999) and immediate 
temperature correction by [Bjorkman & Wood[ ( [200l| ). 

In the second step the level populations are calculated. As dis¬ 
cussed in Appendix [B.2.1[ this is the main problem of line ra¬ 
diative transfer. In the case of circumstellar disks, the problem 
can be significantly simplified by using adequate approximation 
methods. MolSD allows the local thermodynamic equilibrium 
(LTE), the free-escape probability (FEP), or the large velocity 
gradient (LVG) method to be applied. It has been shown that 
LTE is a reliable assumption for lower transitions of abundant 
molecules, but it is recommended that LVG be used, especially 
for high transitions of complex molecules with low abundances 
( [Pavlyuchenkov et al.|2007[ ). 

The last step is to generate velocity spectra or velocity-channel 
maps using a new efficient ray-tracing algorithm. Here, we have 
to consider that the source function is not constant inside a single 
cell, which is due to the underlying velocity field. 

As our code solves the radiative transfer equation including both 
gas and dust, we have to take the optical properties of both into 
account. For the intensity integration, we implemented an em¬ 
bedded Runge-Kutta-Fehlberg solver of order 4(5) with auto¬ 
matic step control. The main advantage is that with providing 
a relative ( 10 “^) and an absolute ( 10 “^^) error limit for the in¬ 
tensity integration, the algorithm calculates the corresponding 
step-width automatically. Conservation of energy is ensured by 
a recursive pixel refinement algorithm based on the underlying 
model grid. 

MolSD is tested and benchmarked against other line and contin¬ 
uum radiative transfer codes. The results are presented in Ap¬ 
pendix]^ along with a detailed derivation of the underlying con¬ 
cepts an^hysical assumptions. The code will soon be publicly 
availabl^under an open-source license. 


3. Model setup 

In this section we investigate the feasibility of observing planet- 
induced structures using molecular lines with ALMA. For 
this purpose we perform 3D hydrodynamic (HD) simulations, 
follow-up line-radiative-transport (LRT) calculations, and simu¬ 
late ALMA observations using the CASA simulator. 


^ https: //github. coin/florianober/Mol3D 
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distance [AU] distance [AU] 

Fig. 1. Example of HD simulations with embedded Jupiter-mass planet (planetmass/starmass-ratio = 10"^). Color-coded is the H 2 number density in 
log scale. Left: Midplane density. The planet has cleared out its orbit, producing a prominent gap. Owing to the perturbations, complex structures 
like spiral arms are also visible. Right: Cut perpendicular to the disk midplane. The gap is clearly visible in the gas distribution. We note that the 
upper layers of the disk have nearly the same density as in the undisturbed case. 


3.1. Three-dimensional hydrodynamic simulations 

In preparation for this study we started with an analytic density 
distribution including a radially symmetric sharp gap completely 
depleted of dust and gas. We find that this approach is insuffi¬ 
cient because it is based on an unrealistically high density con¬ 
trast between the gap and its surroundings, resulting in a gap 
detection rate of almost 100% for the planet-disk configurations 
considered in our study (see Sect. |3.2| ). Instead, justified as¬ 
sumptions about the ratio between the density inside and outside 
the gap, the gap width, and the depth of the gap are needed. 

To reach this goal, we make use of full 3D hydrodynamic sim¬ 
ulations obtained with the well-known finite-volume fiuid dy¬ 
namics code PLUTO ( [Mignone et al.|20()7] ). The code uses the 
HLLC0 solver for the hydrodynamics, a second-order Runge- 
Kutta time integration and second-order variation diminishing 
(TVD) scheme. The setup of the simulations was originally 
presented in [Uribe et aL] (|2011|) and the models used in this 
work are shown in jRuge et al.| ( |2013j ). From their entire set of 
planet-disk configurations we consider two scenarios. The first 
model includes an embedded Jupiter-mass planet (fixed mass ra¬ 
tio Mpianet/^star = 0.001) on a fixed circular orbit. The second 
disk with the same initial parameter setting does not contain a 
planet, and thus serves as an undisturbed smooth reference disk 
for comparison. 

We use a spherical grid where 6 is the angle from the vertical axis 
and 0 is the azimuthal angle. To save computation time, the ver¬ 
tical extent of our PLUTO simulations is chosen to cover only 
6 G [71/2 - 0.3, 71/2 -rO.3], covering four pressure scale heights 
above the disk midplane. 

As initial conditions, it is assumed that the disk gas has a sub- 
Keplerian rotation profile, and that the density profile po oc r~^ 
and the speed of sound Cg = co(rsin6)~^ follow a power-law 
with exponents a = 3/2 and b = 0.5, respectively. The initial 
density distribution is given by 


code 

6 [rad] 

0 [rad] 

K 



PLUTO 

[7r/2-0.3,7r/2+0.3] 

[0,27t] 

256 

128 

256 

Mol3D 

[0,7r] 

[0,27r] 

100 

66 

128 


Table 1. Comparison of grid parameter space of both codes. For the 
original (unsealed) model the radial coordinate ranges from 2 to 9 AU 
for the RT code and from 1 to 10 AU for the PLUTO code because of 
necessary buffer zones. The hydrodynamic simulation results have been 
reduced in order to be usable with Mol3D. 


The disk is described by a locally isothermal equation of state 
P = c^p. The ratio of the pressure scale height h to the radial 
coordinate of the disk is taken to be a constant such that h = 
H/(rsm0) = 0.07. A softened gravitational potential is used for 
the planet. 


(I)(r, 0) = - 


GM 

V|r-rpp + e2’ 


( 2 ) 


where e is the softening parameter which is defined as a fraction 
of the Hill radius e = P rp{Mpl2)^l^ with I = 0.3. 

The original disk configuration has an inner and outer edge of 2 
AU and 9 AU, respectively. The planet is located at = 5 AU. 
For the hydrodynamics simulations, buffer zones are included to 
avoid boundary effects, so the complete radial domain is from 1- 
10 AU. For the follow-up radiative transfer, we need to convert 
the resulting density distribution. We average every two cells in 
each direction to save memory. Thus, we are able to reduce the 
required number of grid cells to For the subsequent radia¬ 
tive transfer (RT) simulations with MoUD we add two cells (one 
above and one below the disk) to fill the entire 6 e [0 , tt] domain. 
The grid parameters for both codes are summarized in Table 


. . ^x-3 (sin^-1 

p(r, 6) = (r sm 6) 2 • exp 


(1) 3.2. Disk model 


Harten-Lax-van Leer-Contact (see 


Toro 2009 1 


The hydrodynamic simulations discussed in Sect. [3.1 [ rely on 
arbitrary units for length and mass. This allows us to introduce 
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linear scaling parameters to scale the disk mass and size to our 
needs for protoplanetary and/or transitional disks. 

Disk size 

For the inner and outer disk radius we chose the linear scaling 
parameter k e [7, ..,25]. Hence, the inner radius of the disk 
is in the range of Rin = k ' 2 AU (i.e., Rin e [14,...,50 AU]) 
and the outer radius in the range of Rout = k • 9 AU (i.e., e 
[63,...,225 AU]. As a result, our disk models mimic different evo¬ 
lutionary states from early protoplanetary disks to more evolved 
transitional disks or circumbinary disks with huge inner cavities. 
The disk density distribution resulting from the simulation of a 
Jupiter-mass planet located at 80 AU in a disk ranging from Rin 
= 32 AU to Rout = 144 AU (corresponding to a scaling parameter 
of k = 16) is presented in Fig. Several specific features in the 
disk density distribution that result from planet-disk interaction 
can be seen. The most prominent one is the distinct gap between 
68 and 92 AU. In the cut through the density distribution per¬ 
pendicular to the disk midplane, the gap appears as a tunnel in 
the protoplanetary disk, as the planet has hardly any influence on 
the disk surface (see Fig. IT] right). Furthermore, the midplane 
density distribution (Fig. [ijleft) shows pronounced spiral arms; 
these features are barely visible in the surface of the disk in this 
case as well. 


Mass range 

We consider a wide range of disk masses (gas-rdust) from 
2.67 -10“^ M© to 2.67 -10“^ M© which approximately covers 
the range of protoplanetary and transitional disk masses derived 
from submm/mm observations (e.g., Andrews et al. [20101 201 1| ). 


Dust properties and stellar parameters 

For the dust we use a mixture of 62.5% astronomical silicate and 
37.5% graphite ( [Weingartner & Draine|200T] ). Furthermore, we 
assume the dust grains to be spherical following a grain size dis¬ 
tribution given in Eq. [3|with an exponent of d = 3.5 (|Dohnanyi| 
[T9691 ): 

dn{a) ~ a~^da. (3) 

We consider dust grains with a minimum radius of amin = 5 nm 
and a maximum radius of a^ax = 250 nm as found in the inter¬ 
stellar medium ( |Mathis et al.|1977| ). The radiation source is the 
central star. For this study we consider a T Tauri and a Herbig Ae 
star with parameters given in Table [^. 


parameter 

T Tauri 

Herbig Ae 

Tstar [K] 

4000 

9500 

Rstar [Rq] 

2 

2.48 

Mstar [M©] 

0.7 

2.5 

Lstar [Lq] 

0.92 

43 


Table 2. Considered pre-main-sequence stars and their properties. 


For the spatial distribution of different molecular species relative 
to each other only weak constraints are available from observa¬ 
tional and numerical studies. For example, observations have 
shown carbon monoxide emission from the cold outer parts of 
protoplanetary disks where CO should be frozen out (<~ 20 K; 
[Dartois et al. |2003[ [Goto et aL]|2012|). A possib le e xplanation 
has been pro posed by | Aikawa & Nomura] ( |2006| ) and [Semenov] 
jet al.| ( |2006| ). They argue, that the CO abundance could be sig¬ 
nificantly higher in the cold (midplane) regions of protoplane¬ 
tary disks owing to vertica l and radial mixing. Motivated by this 
idea, [Hersant et al.j ( |2009| ) computed the chemical evolution of 
circumstellar disks. They included grain surface reactions with 
and without turbulent mixing and CO photodesorption. How¬ 
ever, the resulting CO column densities are still not consistent 
with observations. As indicated by jCleeves et al.j ( |2011| ), the 
disk midplane in the gap region can be irradiated by the stellar 
radiation and X-rays. Thus, CO would not be frozen out onto 
grains, which would result in better gap detection rate thanks to 
higher contrast. 

Because the debate on relative molecular abundances in cir¬ 
cumstellar disks is still open, we decided to use fixed abun¬ 
dances relative to molecular hydrogen. In particular, we fol¬ 
low |P^Iyuchenkovet^ ([20^ and assume their reported abun¬ 
dances for all molecules considered in this study (see Table [^. 
These values are obtained by using a gas-grain chemical model 
with surface reactions ( [Semenov et al.|[2Q04| [2Q05[ ), which is 


based on th e UMIST 95 database of gas-phase reactions ( [Millar 
et al.[[1997[ ). This approach keeps the model simple and hence 


reduces the complexity of our analysis. 


Line radiative transfer setup 

The level populations are calculated by using the LVG method 
(see Appendix [B.2.1[) which has proven its reliability for proto¬ 
planetary disks ( [Pavlyuchenkov et al.[2007[ ). Using these results, 
we produce synthetic velocity-channel maps and spectra for ev¬ 
ery model configuration. We assume a pure Keplerian rotation 
profile. In this study we fix the turbulent velocity Vturb = 100 m/s 
for all of our models, following the results of earlier molecule 
observations (e.g., IPietu et al.|2007l|Chapillon et al.|2012b| ). We 
choose a fixed bandwidth of 50 m/s and calculate 141 velocity- 
channel maps for each of the models with 10° inclination and 71 
for every face-on orientated model in our sample. An overview 
of the considered parameter space of our model is compiled in 
Tablegl 


4. Results 

4.1. Ideal line maps 

In this section we present ideal velocity-channel maps, analyze 
the appearance of planet-induced disk perturbations, and investi¬ 
gate the feasibility of tracing these structures with ALMA (Sect. 
[4.2[ ). In the following we use the term “depth of the gap” as the 
flux density contrast in the velocity-channel maps caused by the 

gap. 


Optical depth effects 


Distribution of the different molecular species 

For the gas-phase, we assume a gas-to-dust-mass ratio of 100:1 
and perfect thermal coupling between gas and dust. 


We begin our discussion based on ideal (3-2) velocity- 

channel maps of our T Tauri disk model. To start with, we 
assume a total disk mass of 2.67 • 10“^ M©, which is at the 
very low end of typical disk masses of these objects (e.g., 
HD 166191, Coka Tau/4; [Kennedy et al.[[MT4[ [Forrest et al. 
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molecule 

transition 

frequency 

[GHz] 

ALMA 

band 

abundance 

N/N(H) 

sensitivity cr 
[mJy] 

PWV 

[mm] 


(1-0) 

115.27 

3 


7.46 

5.2 


(2-1) 

230.54 

6 


3.14 

1.3 

12 c 16 o 

(3-2) 

345.80 

7 

110-4 

4.12 

0.66 


(4-3) 

461.04 

8 


11.66 

0.47 


(6-5) 

691.47 

9 


38.02 

0.47 


(1-0) 

109.78 

3 


3.68 

2.7 


(2-1) 

219.56 

6 


3.11 

1.3 

12cl8o 

(3-2) 

6.2 

7 

2-10-’ 

8.6 

0.66 


(4-3) 

439.09 

8 


61.96 

0.47 


(6-5) 

658.55 

9 


38.66 

0.47 


(1-0) 

89.19 

3 


3.87 

5.2 


(3-2) 

267.56 

6 


3.04 

0.91 

HCO+ 

(4-3) 

356.73 

7 

MO-^ 

4.63 

0.66 


(5-4) 

445.90 

8 


1087.7 

0.47 


(7-6) 

624.21 

9 


137.81 

0.47 


(8-7) 

713.34 

9 


112.10 

0.47 


(1-0) 

88.63 

3 


3.89 

5.2 


(3-2) 

265.89 

6 


2.97 

0.91 

HCN 

(4-3) 

354.51 

7 

M0-® 

4.48 

0.66 


(5-4) 

443.12 

8 


46.78 

0.47 


(7-6) 

620.30 

9 


4.89-10^ 

0.47 


(8-7) 

708.88 

9 


59.91 

0.47 


(2-1) 

97.98 

3 


3.38 

2.7 


(3-2) 

146.97 

4 


3.17 

1.8 


(5-4) 

244.94 

6 


2.9 

0.91 


(6-5) 

293.91 

7 


3.82 

0.91 

CS 

(7-6) 

342.88 

7 

M0-* 

4.12 

0.66 


(8-7) 

391.85 

8 


9.77 

0.47 


(9-8) 

440.80 

8 


27.2 

0.47 


(10-9) 

489.75 

8 


32.11 

0.47 


(13-12) 

636.53 

9 


35.96 

0.47 


(14-13) 

685.44 

9 


33.13 

0.47 


Table 3. Properties of the different molecules considered in this study. We note that the sensitivity is calculated with the ALMA sensitivity 
calculator (version: April 2015) for an observation of three hours considering 50 antennas in Jwa/-polarization mode. The sensitivity strongly 
depends on the atmospheric transmission, which explains the high values for some transitions. We use the recommended values for the precipitable 
water vapor (PWV) column density for each transition, which are also included in the sensitivity calculator. 


parameter 

value(s) 

Rin [AU] 

14-50 

Rou, [AU] 

63 - 225 

[Mq] 

2.67 • 

Pdust [g/cm^] 

2.5 

3^min 

0.5 

^max [A^I^] 

250 

inclination [°] 

0, 10 

Vturb [m/s] 

100 

distance [pc] 

140 

observing time [h] 

3 

baseline [km] 

0.7 - 16.3 

bandwidth Av [m/s] 

50 


Table 4. Parameter space of the disk model. 


|2004[ ID’Alessio et al.||2005| ). The inner rim of the chosen disk 
has a radius of 32 AU while the outer radius amounts to 144 AU. 


The planet is located at ~80 AU and the planet-induced gap has 
a width of ~32 AU. The disk is seen face-on. In this case, the 
line shows a symmetric gap throughout the different velocity 
channels, because the line width is given only by thermal and 
microturbulent velocities (see Fig. see also Appendix |B.2| ). 
We find that the contrast between the gap and its edges is not 
constant throughout the individual velocity channels, but varies 
with the optical thickness of the chosen line. In the optically 
thin case, the gap depth increases as it approaches the line center 
owing to the increasing line flux. In our case, the (3-2) 

transition is optically tick at its center even for this lowest disk 
mass of our sample. Thus, the r = 1 surface seen in the line 
wings belongs to parts of the disk closer to the disk midplane 
than the disk layers traced at the line center (i.e., at 0 m/s). 
Therefore, for this specific configuration we observe the gap 
best at velocities + 300 m/s, while the main contribution at the 
line center comes from the upper layers of the disk. As these 
upper layers are less disturbed by the planet-disk interaction, 
the gap is less pronounced in the center of an optically thin line, 
i.e., the nearly undisturbed upper disk layers mask the gap. In 
analogy, further signatures of the planet-disk interaction like the 
pronounced spiral arm in the outer disk can be identified best in 
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Fig. 2. Ideal (3-2) velocity-channel map of a T Tauri disk model with an embedded Jupiter-mass planet with an orbital distance of 80 

AU. The disk mass amounts to 2.67 • 10"^ M©, which is the lowest disk mass in our sample, and the outer radius is 144 AU. The flux density in 
Jy/mas^ is color coded for each individual channel. The disk is seen face-on (0°). The gap is clearly visible, particularly in the line wings. Even 
more complex features like the huge spiral arm (e.g., at 300 m/s) can be recognized. The total width of the line is only given by thermal and 
microturbulent velocities (Section [R^. 


the line wings (here, in the + 400 m/s channel maps). 

If we increase the optical depth by increasing the disk mass by 
two orders of magnitude (i.e., 2.67 • 10“^ Mq), which represents 
a typical disk mass of protoplanetary disks, we find that the 
gap vanishes nearly completely in the velocity-channel maps 
(Fig. [g. As only the envelope above the gap can be observed, 
the disk appears to be undisturbed without any clear signs 
of planet-disk interaction. This is the main problem when 
observing For disk masses corresponding to typical 

protoplanetary disks masses we find that the molecule is 

not an adequate tracer for gaps. In this case, other less abundant 
molecules/isotopologues that allow deeper disk layers to be 
observed are better suited. We illustrate this on the basis of the 
(3-2) transition (Fig. |^. Because of the low abundance 
of this carbon monoxide isotopologue, this line is less optically 
thick and thus traces deeper disk layers. However, its excitation 
temperatures of the lower energy levels are very similar to those 
of the molecule. In contrast to the optically thick 

line of the same model, the gap is clearly detectable. 

The most massive disks (2.67 -10“^ M©) in our sample do not 
show any signs of a gap in any of the considered transitions (see 
Table 1^. For the reasons outlined above, for a specific molecule 
observation with the aim to reveal gaps, we recommend estimat¬ 
ing the optical thickness of the disk (for example from its mass 
derived through continuum observations) in order to choose the 
appropriate molecule and (optically thin) transition. 


Disk inclination 

In contrast to the ideal face-on models, the velocity-channel 
maps of a disk with a small inclination of 10° is shown in Fig. 

Negative velocities (top left corner) belong to parts of the 
disk moving toward the observer and areas with positive (bot¬ 
tom right) velocities are moving away from the observer. In this 
figure the main problem for detecting gaps in inclined disks, and 
thus the majority of potential targets, becomes obvious. Only 
fractions of the gap can be seen, but it does not appear as a sym¬ 
metric ring as the gas covers a broad range of velocities (» Vturb) 
along the gap throughout the individual velocity channels. This 
effect is illustrated in Fig.[^ The underlying disk model (Fig.[^ 
left) has an outer radius of 117 AU and the planet is located 
at 65 AU. Owing to the disk inclination of 10°, the line is sig¬ 
nificantly Doppler shifted. The velocity-channel map shows a 
cavity in the orbital region of the planet. However, this not an 
irrevocable sign of a gap. Owing to the projected velocity in the 
observer’s direction, the undisturbed model (see Fig. right) 
also reveals a prominent local minimum in the same disk region. 
Thus, the velocity patterns effectively hide the gap and hinder 
unambiguous identification, especially for optically thick lines. 
A possible solution is to use velocity integrated intensity maps 
for gap detection instead (see Fig. [7]). For low disk inclinations, 
these images do not deviate much from face-on inclined disks, 
i.e., for an inclination of 10° as considered in this study the flux 
deviation compared to the disk model seen face-on is less than 
5%. The velocity integration has the advantage that the patterns 
disappear and the gap can be identified easily. However, the ve- 
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Fig. 3. Ideal (3-2) velocity-channel map of a T Tauri disk model with an embedded Jupiter-mass planet (at 80 AU). The disk mass 

amounts to 2.67 • 10"^ M© and the outer radius is 144 AU. The flux density in Jy/as^ is color coded for each individual channel. The disk is seen 
face-on (0°). The gap is barely visible. The line is optically thick owing to the material above the gap. However, the flux density is comparable to 
the model with 100 times lower mass (Fig. |^. 



Fig. 6. HCO^ (4-3) velocity-channel map at 150 m/s. The disk has a 
slight inclination of 10°. Disturbed (left) and smooth disks reveal simi¬ 
lar results. Both models show indications of a gap. Thus, unambiguous 
gap detection requires detailed analysis of the origin of those structures. 
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Fig. 7. (3-2) velocity integrated intensity map of the channel 

map shown in Fig.[^ Owing to the integration, the patterns due to the 
underlying velocity held have vanished. 


locity dependency of the gap depth (see previous section) is lost. 
Consequently, the velocity-integrated maps are dominated by the 
line flux at the line center, which traces only the upper disk lay¬ 
ers in the case of optically thick lines. 


Central star 

We now discuss the impact of using the Herbig Ae star, which is 
~ 45 times more luminous, as the central star. 

As the disk temperature is higher than the T Tauri case, the 
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Fig. 4. Ideal (3-2) velocity-channel map of a T Tauri disk model with an embedded Jupiter-mass planet (at 80 AU). The disk mass 

amounts to 2.67 • 10"^ M© and the outer radius is 144 AU. The flux density in Jy/as^ is color coded for each individual channel. The disk is seen 
face-on. Mostly because of the low abundance of this CO isotopolo^e, this line traces deeper parts of the disk. Thus, the gap is clearly detectable 
in contrast to the optically thick line of this model (cf. Fig. |^. 


molecules are exited to higher energy states. The line radiation 
of the outer parts of the disk increases because a larger fraction 
of molecules is sufficiently exited to the relevant energy levels 
due to the warmer disk. In contrast to the situation in the outer 
disk regions, the higher temperature in the innermost region re¬ 
sults in a depopulation of the lower energy states. Consequently, 
the low line transitions we consider in this study produce less 
flux in this region. 

Second, the mass of the Herbig Ae star is about 3.5 times higher, 
resulting in a higher Keplerian velocity at a given radial distance 
from the central star. Thus, for the inclined disks the linewidth 
of the Doppler broadened line is much higher than in the T Tauri 
star case. Consequently, at a certain velocity channel (band- 
widths of 50 m/s in this study) the parts of the disk with the 
appropriate velocities to contribute to this channel appear pre¬ 
cise and well defined (see Fig. [^. Thus, assuming the same 
disks for both PMS stars, individual disk regions can be better 
studied in the case of the more massive central star, which in¬ 
creases the chances of gap detection. This effect is a direct result 
of the dependence of the linewidth on the kinetic temperature 
(cf. Eq.[^. 


4.1.1. Detectability of gaps in the ideal case 

We describe our approach to identifying a gap and to deriving its 
depth in the velocity-channel maps. 

We follow the procedure outlined by |Ruge et ^ ( |20I3| ). We 
derive the contrast between the inner/outer edge of the gap and 
the gap itself by using radial cuts throughout the simulated maps 
(see Fig. [^. To reduce the probability of false detection (e.g., 
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Fig. 9. Illustration of the method of detecting gaps. The blue line 
indicates a radial cut through the HCO^ (4-3) velocity-channel map 
(Fig. [14} at the central frequency (v = 0.0 m/s). For the simulated 
ALMA observations, the gap is considered as being detected if the depth 
of the gap is deeper than three times the noise level cr of the selected 
transition. In this particular case, the gap is clearly detectable, resulting 
in a gap contrast of ~ 7cr. 


local minimum in the brightness distribution due to the inclina¬ 
tion effect discussed in Sect. |4.1| ), 64 individual radial cuts are 
considered in the analysis. Only if a gap is indicated along more 
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Fig. 5. Ideal (3-2) velocity-channel map of a T Tauri disk model with an embedded Jupiter-mass planet. The flux density in Jy/as^ is 

color coded for each individual channel. The inclination of the disk is 10° and the disk has a total mass of 2.67 • 10“^ M©. The gap is visible but 
does not appear as a symmetric ring like in the continuum case because the parts of the disk contributing in the individual channel maps change 
due to Doppler shift. 



than 50% of successive radial cuts is it considered as detected. 
We successfully tested our algorithm with smooth undisturbed 
reference hydrodynamic models to verify that no substructure is 
misleadingly considered as a gap. However, because of the very 
general approach due to the huge parameter space, we emphasize 
that for specific observations better results may be achieved by 
adjusting the algorithm to the specific observation (e.g., a well- 
adapted compensation of background noise/sparse uv-coverage, 
which leads to asymmetric structures). 


Overview charts 

In this section, the feasibility of detecting gaps through obser¬ 
vations of selected molecular lines for all considered disk con¬ 
figurations (see Table is summarized. For this purpose, we 
derive the gap depth from the ideal velocity-channel maps for all 
five molecules and 32 assumed transitions using the algorithm 
described above. For every transition we derive the expected 
ALMA background noise of the reconstructed velocity-channel 
maps (cr) using the ALMA sensitivity calculator (see Sect.|4.2.I| 
Tablej^. Owing to the atmospheric transmission, cr mostly rises 
with frequency resulting in a worse signal-to-noise ratio. To take 
care of this effect, we derive the gap depth and weight it by the 
expected sensitivity of ALMA at every line transition considered 
in this study. Thus, we obtain a criterion of how well individual 
transitions allow gap detection with ALMA. These results allow 
those transitions to be chosen that are best suited for the simu¬ 
lated ALMA observations in Sect. 14.21 

Figure [T^ shows the results for an exemplary disk with an outer 


radius of 144 AU and a disk mass of 2.67 • 10 ^ M© around the 
Herbig Ae star. The gap depth measured in every velocity chan¬ 
nel for the (3-2), (3-2), and CS (7-6) transition 

is shown. All transitions are optically thick at the line center. 
While the (3-2) and CS (7-6) transitions fail to detect the 

gap at the line center, the optically thin (3-2) transition 

allows for gap detection even at the central frequency. Neverthe¬ 
less, the highest gap depth is achieved by the CS (7-6) transition 
at the line wings (at +350 m/s). 

To summarize our findings for all molecules, we search for the 
transition that provides the highest gap depth for given disk mass 
and disk size and present the result in one overview map. The 
resulting charts for both considered PMSs are shown in Fig. 

The full overview charts showing the individual gap depth of all 
32 transitions considered in this study can be found in Appendix 

El 

As already indicated in Sect. |4.I[ the most promising molecule 
for identifying gaps in low-mass disks (-lO""^ - 10“^ M©) is 
The gap can be detected either at the line center for the 
lowest disk masses or at the line wings for the intermediate-mass 
disks. In most configurations the fiux density and consequently 
the gap depth is higher even in the line wings compared to the 
gap depth which can be achieved with other molecules. 

For the disks with intermediate masses (~I0“^ - 10“^ M©) of 
our sample gap detection becomes very challenging. In these 
cases ^^C^^O is optically thick for almost all considered disk 
configurations even at the line wings. Thus, transitions of less 
abundant molecules like HCO^, HCN, or CS are better suited 
for gap detection. However, because of the low abundance of 
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Fig. 8. Ideal (3-2) velocity-channel map of a Herbig Ac disk model with an embedded Jupiter-mass planet. The disk mass amounts to 

2.67 • 10"^ Mq and the outer radius is 144 AU. The flux density in Jy/as^ is color coded for each individual channel. The inclination of the disk 
amounts to 10°. 
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Fig. 11. Overview plot of the maximum gap depth in the ideal velocity-channel maps when comparing all transitions considered in this work. 
The central radiation source is a T Tauri (left) or a Herbig Ac star (right). For the most massive disks in our sample. All considered transitions 
are optically thick and only trace the surface of the disks. Thus, for the smallest disks the gap is detectable outside the line in the dust continuum 
(marked with Dust). 
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Velocity [m/s] 


Fig. 10. Depth of a gap for three different molecules/transitions over 
velocity and weighted by the ALMA noise level cr. The disk has an 
outer radius of 144 AU and the disk mass amounts to 2.67 • 10"^ M© 
and the disk heated by the Herbig Ae star. For this disk model, all con¬ 
sidered transitions are optically thick at the line center and the detection 
of the gap is only feasible with due to its low abundance. Nev¬ 

ertheless, at the line wings other molecules transitions allow for gap 
detection too, and in this particular case the gap depth has a maximum 
value when using the CS (7-3) transition. 


these molecules the net flux is very low. Consequently, the gap 
depth that relates the flux density to the background noise is too 
low for gap detection in many cases. The most massive 
Mo) disks in our sample do not show any reliable sign of a gap 
for any molecule transition considered. We find that the gaps of 
these disks can be traced best outside the line in the dust contin¬ 
uum, even with the low bandwidth considered in this study. 
When comparing the molecular line emission from the disk 
around both PMSs, we find very similar results in terms of the 
probability of gap detection. The effects due to their different 
stellar mass and luminosity (see Sect. |4.1| for a discussion) lead 
to slightly different molecules, which offer the best probability 
of detecting a gap (e.g., because of different excitation states 
allows for gap detection even for the disks with a mass 
of ~10-4 Mo). 

Therefore, our results demonstrate the feasibility of gap detec¬ 
tion using molecular lines, but it is only efficient for low-mass 
protoplanetary disks or disks with huge inner cavities (e.g., tran¬ 
sitional disks or circumbinary disks). 


4.2. Simulated ALMA observations 


4.2.1. ALMA simulation setup 


We use the Common Astronomy Software Application pack¬ 
age (CASA 4.2, |Petry~ fe CASA Development Team|[2012| ) 
to simulate observations based on our ideal velocity-channel 
maps. We consider a total of six medium and extended con¬ 
figurations (baselines from 0.7 to 16 km) using 50 antennas, 
a total observing time of 3 hours, a target distance of 140 pc, 
and the sky-coordinates of the “Butterfly Star” (a = 04^33”^ 17^ 
S = -r22°53'20”, |Wolf et al.|2003l|Grafe et al.|2013| ) which is 


representative of young stellar objects in the Taurus-Auriga star¬ 
forming region. The precipitable water vapor (PWV) in the at¬ 
mosphere has a major influence. Therefore, we also include ther¬ 
mal noise in our simulations and use the recommended column 


densities for the line transitions we consider in our sample (see 
Table [^. First, we employ the simobserve task to generate real¬ 
istic visibility measurements. The final simulated observations 
are obtained using the simanalyze tas k which uses th e CLEAN- 
algorithm for image reconstruction ( |Hogbom||1974| ). The ex¬ 
pected noise levels cr are estimated with the ALMA sensitivity 
calculatoij^ because the measured noise in the simulated images 
strongly depend on the quality of the deconvolutioij^ The prop¬ 
erties of the considered molecules and the expected sensitivity 
are summarized in Table We consider the ALMA polariza¬ 
tion configuration in dual-mode. With respect to the results from 
Sect. |4. 1.1 1 we select those molecules that offer the highest prob¬ 
ability of gap detection, i.e., (3-2), (3-2), CS (7- 

6), HCO^ (3-2), andHCN (3-2). In addition, their corresponding 
transitions around 345 GHz (ALMA band 7, see Table [^, of¬ 
fer a reliable comparability of the different molecules and yield 
a good compromise between spatial resolution (higher frequen¬ 
cies) and better atmospheric transmission (lower frequencies^ 
The (4-3) transition is added to provide a second 

transition offering higher spatial resolution at the cost of sensi¬ 
tivity. 


4.2.2. Simulated line observations 

In Fig[T^ and Fig[^ a face-on and an inclined disk model are 
shown, respectively, as they could be observed with ALMA us¬ 
ing configuration 14 (~ 1.6 km baseline) within three hours (cor¬ 
responding ideal maps are shown in Figs, [^an d0 respectively). 
More compact baselines result in higher flux densities and con¬ 
sequently higher signal-to-noise ratios, while larger baselines al¬ 
low gaps of compact disks to be traced at the cost of decreasing 
signal-to-noise ratio. We find for the most extended disk con¬ 
figurations (Rin ~ 45 AU, Rout ~ 200 AU) in our model setup, 
maximum baselines of about 1 km are sufficient to identify the 
gap in these disks, but fail for the smaller disks. 

For the compact disks in our parameter space, higher spatial res¬ 
olution and thus longer baselines are required. In our sample the 
smallest disks have an inner/outer radius of 14 AU/63 AU, re¬ 
spectively. We find that unambiguous gap detection is feasible 
using intermediate baselines of about 1-3 km. Longer baselines 
are not required and only result in sparser uv-coverage and thus a 
lower S/N of the reconstructed maps. Furthermore, many of the 
smaller planetary signatures (inner and outer spiral arms and disk 
inhomogenities) are eliminated/smoothed out because of beam 
convolution and worse uv-coverage. However, for this particu¬ 
lar model, and for most of the expanded models in our sample, 
the gap can be identified in both inclination cases using ALMA 
configuration 14. For the same disk, the resulting HCO^ (4-3) 
velociW-channel map for ALMA configuration 14 is shown in 
Fig. The HCO^ molecule is less abundant than and 

consequently the line is less optically thick. Therefore, regions 
closer to the disk midplane can be observed and the gap appears 
more pronounced. In contrast, the line flux itself is about 35% 
weaker than the simulated (3-2) line of the same model. 


^ see https://almascience.eso.org/proposing/ 

sensitivity-calculator version April 2015 

see http: //casaguides.nrao. edu/index.php?title=Guide_ 
To_Simulating_ALMA_Data 

see https://almascience.eso.org/about-alma/weather/ 
atmosphere-model 
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Fig. 12. (3-2) velocity-channel map of a model with an embedded Jupiter-mass planet as it could be observed using ALMA configuration 

14 (1.6 km baseline). The inclination of the disk is 0° (face-on). The outer disk radius amounts to 144 AU and the disk mass to 2.67 • 10"^ M©. 
The flux density in Jy/beam is color coded for each individual channel. In contrast to the ideal velocity-channel maps, complex structures are 
eliminated owing to beam convolution and noise. However, the gap can still be distinguished at the line slopes (+ 300 m/s). 
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Fig. 13. (3-2) velocity-channel map of the same model shown in Fig. 

is color coded for each individual channel. 
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but with a disk inclination of 10°. The fiux density in Jy/beam 
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Fig. 14. HCO^ (4-3) velocity-channel map of a model with an embedded Jupiter-mass planet as it could be observed using ALMA configuration 
14 (1.6 km baseline). The flux density in Jy/beam is color coded for each individual channel. The inclination of the disk amounts to 0° (face-on). 
Compared to the same model but simulated (3-2) observation, the line flux is about 35% lower owing to the distinct lower abundance of 

HCO^. However, because of the low abundance the gap is clearly detectable, even at the central wavelength. 


4.2.3. Gap detection results 


Our main intention is to investigate the feasibility of tracing gaps 
in molecular lines with ALMA. To allow this for a large range of 
disk models, we define a general criterion for a quantitative mea¬ 
sure of gap detection in the simulated ALMA maps. W e apply 
the same algorithm as in the case of ideal maps (see Sect. |4.1.r] ). 
The resulting gap depth is given in units of cr to clarify the mean¬ 
ing of the gap detection. Only if the depth of the gap is greater 
than three times the noise level cr do we consider the gap as de¬ 
tected. 

First, we demonstrate that the optical depth effect (see Sect. |4.1| 
for a discussion of the ideal case) has direct observational im¬ 
pact. As shown in the ideal case, the choice of the molecule for 
a given disk mass/size is most important for unambiguous gap 
detection. In Fi g.p~5] we present the gap depth for the same tran¬ 
sitions as in Fig.|10[ but now for a simulated ALMA observation 
utilizing configuration 14 (~ 1.6 km baseline). Both results are 
closely comparable as the gap depth in the ideal case shows the 
same trend as the simulated ALMA observations, i.e., ^^C^^O(3- 
2) cannot be used as a tracer for gaps, and (3-2) allows 

for gap detection even at the line center, but the CS (7-6) transi¬ 
tion traces the gap best and results in a gap depth of about 11 cr 
at + 350 m/s. 

Second, our results are summarized in the overview charts for 
the Herbig Ae (left) and T Tauri (right) star (Fig. [^. Compact 
ALMA configurations (9, 11) result in high signal-to-noise ratios 
(S/N), but lack of spatial resolution. As a consequence, utilizing 
the most compact ALMA configuration nine gaps are only iden¬ 
tifiable for the most extended models of our sample (e.g., disk 



Velocity [m/s] 


Fig. 15. Measured gap depth in units of cr over velocity. Gap detec¬ 
tion with the optical thick (3-2) transition (blue) is feasible, but 

only at the line edges (e.g., + 700 m/s). In contrast, the less abundant 
(3-2) transition (green) allows for gap detection at the line cen¬ 
ter. However, the CS (7-6) gives the most significant gap depth of about 
17cr at + 350 m/s. 


models with Rout = 255 AU). If we consider longer baselines 
(higher spatial resolution) in the range of -1.6 km - 2.3 km we 
are able to identify gaps for more compact models with a cr > 3. 
Intermediate baselines (e.g., configuration 11-14, -1.0-1.6 km) 
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Fig. 16. Overview maps of the detectability of gaps around Herbig Ae (left) and T Tauri (right) stars. Each column represents one ALMA 
configuration and each row represents one molecule transition. The x-axis of the inner plot gives the total mass of the disk, and on the y-axis the 
scaling factor k is linearly connected to the spatial dimension of the disks. The depth of the gap is color coded in units of cr. We only consider 
a gap as detected when the depth is at least three times the background noise cr. For both stars, gap detection is feasible for a wide range of 
models in our setup within three hours of observation. The difference between the two stars is mostly due to the much higher luminosity. Similar 
to the results in the ideal case, gap detection is feasible but is restricted to low-mass disks (M^isk < 10"^). Gap detection using the most extent 
configuration (max. baseline ~16 km) is not feasible within three hours of observation, which is mostly due to the small bandwidth of 50 m/s and 
is not shown in this overview plot. 


potentially allow gaps to be traced for the wide range of models 
considered in this study. 

The comparison between the molecules considered shows that 
CS, HCN, and HCO^ allow the detection of gaps even for disks 
with intermediate masses (~10“^ M©), while tends to 

fail. As already discussed in Sect. |3.2| owing to the optical thick 
12 c16o line only the surface of the disks can be observed, which 
is barely perturbed by the planet-disk interaction. In contrast, for 
lines of molecules with lower abundances the surface is still op¬ 
tically thin even for more massive disks. We find that especially 
CS and HCN are proper candidates to serve as a tracer for gap 
detection in the mass range from ~10“^ M© to 10“"^ M©. 

For more massive disks (Mdisk > 10“^ M©) the ^^C^^O isotopo- 
logue ^^C^^O is an adequate candidate. However, owing to the 
resulting low line flux gaps these massive disks should be ob¬ 
served in the continuum because in the submm wavelength range 
provided by ALMA the di sk is still optically thin in the gap re¬ 
gion (see |Ruge et al.|2013] ). 

Finally, we And that because of the assumed observation time 
of 3 h, long-baseline configurations (20, 28) are not suited to 
identifying gaps because of the low intrinsic flux and inadequate 
uv-coverage. 


5- Summary 

Our intention was to explore how planet-induced structures, 
in particular planetary gaps, appear in high-angular resolution 
molecular line maps, and subsequently if ALMA would allow 
the detection of these structures. For this reason, we devel¬ 
oped the line radiative transfer code called MoUD. It allows 
the calculation of the level-populations with several N-LTE ap¬ 
proximate methods and the production of ideal velocity-channel 
maps. First, we tested the new code against other established 
(molecular line) RT codes. 

Based on 3D hydrodynamic simulations carried out with the 
PLUTO code, we simulated synthetic velocity-channel maps for 
five molecules HCO+, HCN, and CS) for a to- 

tal of 32 rotational transitions. We summarize our findings: 

- The patterns of the global (Keplerian) velocity field in the 
corresponding channel maps of inclined disks mask planet- 
induced structures. Gap identification can be ambiguous in 
these cases. 

- Molecules that tend to become optically thick (like 

only trace the upper layers of the disk, which are rarely af¬ 
fected by the planet-disk interaction. In consequence, gap 
detection for massive disks (> 10“^ M©) is barely feasible 
with molecular lines. 
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- For lower mass disks, depending on the molecule and its op¬ 
tical properties (abundances, excitation state) the gap can be 
observed for a broad range of considered disk configurations, 
not necessarily at the central wavelength of this transition 
but at the line wings. Because of the higher optical depth of 
the line at its center, the r = 1 surface traces only the upper 
unperturbed disk layers, while the r = I surface of the line 
wings traces deeper and thus more strongly perturbed layers. 

- With ALMA it is feasible to identify gaps using molecular 
lines for several disk configurations, if one respects the ef¬ 
fects we outlined above (i.e., by choosing an appropriate op¬ 
tically thin molecule and corresponding transition). In par¬ 
ticular, for low-mass disks (< 10“^ M©) the feasibility of 
identifying gaps using molecular lines is very promising. For 
this mass range molecular line observations fill the discovery 
space for gaps where continuum observations are supposed 
to fail (see Fig. 8 in the complementary study of|Ruge et al.| 
[20T3l ). 

- Our simulations are based on an assumed distance of 140 pc 
and a fixed sky coordinate (max. elevation ~ 45°) to give 
reliable results for the Taurus star-forming region. Thus, for 
many disks located closer or at more favorable positions on 
the sky the feasibility of detecting gaps is expected to be even 
higher. 


We emphasize, that the observation of molecules with different 
abundances allow for unambiguous identification of gaps and the 
physical and chemical composition of protoplanetary disks. Es¬ 
pecially for low-mass disks where gap detection in the contin¬ 
uum is difficult owing to the low intrinsic flux, molecular line 
observations are a promising alternative with which to identify 
gaps. 

Finally, we note that our study is based on HD simulations. In 
addition, |Uribe et ah] ( |2011| ) have remarked that magnetohydro¬ 
dynamics (MHD) simulations result in more extended gaps ow¬ 
ing to the presence of the weak magnetic field. Including intrin¬ 
sic stellar radiation, the form of the gaps also changes because 
the gap walls can be directly radiated by the host star, result¬ 
ing in deeper gaps. This effect is even higher for the outer wall 
( |Jang-Condell & Turner|20I2|). Using smoothed par ticle hydro- 
dynamic (SPH) simulations, [Gonzalez et al.| pQ12| ) found that 
gaps are deeper for larger dust grains and bigger planets. In ad¬ 
dition, we only considered a Jupiter-mass planet in this study. 
More massive planets would enlarge the gap width owing to 
the higher gravitational potential. If these effects are present it 
would therefore be even easier to observe gaps with ALMA than 
proposed within this study. 
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Appendix A: Gap depth - complete overview charts 

The full overview charts for all molecules considered in our 
stud y are summarized in Fig. |A.1| for the T Tauri case and in 
Fig. |A.l| for the Herbig Ae case. The derived cr-weighted gap 
depth is color coded for each considered disk configuration. 


Appendix B: Mol3D 

In this section we present a detailed derivation of the underlying 
methods, physics, and assumptions used by our new N-LTE line 
and dust continuum radiative transfer code Mol3D. 

Appendix B.1: Radiative transfer equations 

We start with the transfer equation in the following form, 

= -ay(s)Iy(s) + jy(s), (B . 1 ) 

with the monochromatic intensity /y, the emission factor jy(s), 
and the total absorption factor ay(s). The intensity from point pi 
to p 2 along its path s can be calculated using the integral form of 
equation [B.l[ 

Iy(P2) = + r Jy(s)e-^’'^^’P^^ds, (B.2) 

Jpi 

where Ty is the optical depth between two points, pi and p 2 , 
along the ray of sight defined as 

^P2 

Ty(pi,P2)= I a'v(s)d^. (B.3) 

Jpi 
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Appendix B.2: The gas-phase 

For the gas, e.g., for a selected transition from lower level i to 
upper level j, the emission (jy) and absorption {ay) coefficients 
have the following form: 


hv 

jijiv) = -^NniAij<Pij{v), 

hv 

aij(v) = —NitijBji - 


(B.4) 

(B.5) 


We follow the ansatz of |Rybicki & Hummer ( 1991| l for the non- 
overlapping multi-level line treatment. For a molecule with N 
levels and a selected transition between lower i and upper j level, 
we consider the spontaneous downward transition rates Aij, the 
Einstein coefficients Bij for stimulated transition, and the colli¬ 
sion rate Ctj. The collision rate can be calculated from the num¬ 
ber density of the collision partner and the downward collision 
rate coefficient yij, which is the Maxwellian average of the col¬ 
lision cross section cr measured in laboratory: 




ij ~ ricoVYij' 


(B.6) 


All molecule data we use in our software package (e.g., frequen¬ 
cies, transition rates Aij, collision rate coefficients) are obtained 
from the Leiden Atomic and Molecular Database (LAMDA) 
( Schoier et al.|2005] ). As a consequence, Mol3D uses the same 
common input format for the molecule data as many other line 
radiative transfer programs a vailable like RADEX ([va n der Tak 
et al.pOOT] ) or URAN(IA) ( [Pavlyuchenkov & Shustov ]2004| ). 
This approach offers the advantage of easily extending the 
code for several common molecules available in the LAMDA 
database (3 atomic and 33 molecular species at the time of writ¬ 
ing). The dominant broadening effects are thermal and turbulent 
broadening. Hence, we assume a Gaussian line profile function, 
but in principle any other profile function can be included as 
well. 


0,y(y) 


VtotV/y 


exp 


2 2 

Vii 


(B.7) 


where Vij is the central frequency of the transition and c the speed 
of light. The parameter Vtot gives the total line broadening factor 
taking into account the assumed microturbulent velocity Vturb and 
the kinetic velocity Vkin due to the kinetic temperature Tkm: 


Vtot 


= V' 


V^ -|- 
'^kin ^ '^turb 


2^rkin 

rrimoX 


+ v; 


turb’ 


(B.8) 


In our code, the microturbulent velocity Vturb can either be set to 
a common value, ~ 0.1 - 0.2 km/s for protoplanetary disks ( jPietu 
et al.|2007} [Hughes et al.pOlT] ), or for example in the case of an 
outcome of (M)HD/MRI simulations, it can be provided as an 
input for every grid cell. 

Finally, to describe the full line radiative transfer problem, a set 
of balance equations is needed: 


^ y rijAij -b {niBij njBji)Jij] ^^{njAji {njBji niBij)!p 


i>j 


KJ 


(B.9) 




■ - 0- 


In fact, Eq. B.9 is not one single equation, but a whole set of 
linear equations, one for every energy level. In our code this lin¬ 
ear matrix equation is solved with a simple Gaussian elimination 
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Fig. A.l. Overview plot of gap depth in the ideal velocity-channel maps for five different molecules. The central radiation source is a TTauri 
star. The cr-weighted gap depth is color coded. 


algorithm. This set has to be solved locally (for every grid cell), 
but it has a global character due to its dependency on the line 
integrated mean intensity Jij defined as 




1 

471 



Iy(^(y)dy. 


(B.IO) 


Equations |B .9 1 and [B . 1 0| are directly coupled to Eq. |B.1[ 

To solve the line radiative transfer problem, one needs to start 
with estimated values for the level populations, solve the ra¬ 
diative transfer equation, calculate the mean intensity, and then 
solve Eq. |B.9| again to create a new set of level populations. 
Then one needs to iterate until the correct level populations have 
been found. 


|et al.| ( [2007| ) have shown for the case of protoplanetary disks, 
for the most common molecules LVG is a very accurate method 
of calculating the level populations and gives a very good ratio 
between computational efficiency and reliability of the results. 


LTE level populations 

In this assumption the level populations follow a Boltzmann dis¬ 
tribution and the excitation temperature of the selected transition 
is assumed to be equal to the kinetic temperature Texc = T'kin- 


«i / hvij 

— = — exp -- 

Hj gj \kTun 


(B.ll) 


Appendix B.2.1: Level populations 

The main problem in radiative line transfer is the calculation of 
the level populations. Equation ( |B.9| ) shows that the level pop¬ 
ulations are directly coupled to the mean intensity and therefore 
to the gas temperature, which is unknown for complex structures 
as in the case of protoplanetary disks. 

We implemented dissimilar algorithms to calculate the level 
populations with three approximate methods, for example, lo¬ 
cal thermodynamical equilibrium (LTE), full escape probabil¬ 
ity (EEP), and large velocity gradient (LVG). As |Pavlyuchenkov| 


Here gi and gj are the statistical weights of the i-j line transition 
and yij the corresponding central frequency. 


LVG level populations 

In the case of protoplanetary disks, the radial velocity gradients 
are usually much larger than the local thermal and microturbu- 
lent velocities ( |WeiB et aL |2005[ [Castro-Carrizo et J^|2007| ). 
This means that photons emitted at a certain disk region can only 
interact and consequently get absorbed locally. In this approach, 
the mean intensity J is approximated from the local source func- 
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Fig. A.2. Overview plot of gap depth in the ideal velocity-channel maps for five different molecules. The central radiation source is a Herbig Ac 
star. The cr-weighted gap depth is color coded. 


tion S in combination with an external radiation field /ext- 


let al.|2007H : 


J = {\ -li)S + phn, 

0<y8< 1. 


(B.12) 

(B.13) 


T(y) = a{v)R 



(B.15) 


The quantity p gives the probability of a photon beeing able to 
escaping the model. We note that in the optically thin case a 
photon can escape the model without interaction. For this case, 
P equals 1 and the LVG method would be equal to the full escape 
probability method. Therefore, the FEP method (also included 
in Mol3D) uses the same formalism as the LVG method, but with 
P always set to 1 (see also [van der Tak et al.|2Q07| ). To calculate 
the local p, our program uses Eq . B.14 introduced by |Mihalas| 
|et al.| ( p778) ); |de Jong et al.| ( |1980| ) for accretion disks in general. 




1 - exp(-T) 


(B.14) 


where r is the effective optical depth of the observed line. 
We note that this formula is well suited in the case of proto- 
circumstellar disks. For other geometries and scenarios other 
approximations might give better results (see, e.g., |de Jong et al.| 
( 1975| ) or IQsterbrock & Ferland] ( |2006| )). At this point an ap¬ 
proximation of the optical line depth is made. As it is connected 
to the local velocity field V, microturbulence Vturb. and gas den¬ 
sity, we use the following expression (see also |Pavlyuchenkov| 
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Appendix B.3: Benchmarks 

We demonstrate the applicability of the code Mol3D. For this 
purpose we compare its results to t hose obtained with the line 
radiative transfer code URAN(IA) ( |Pavlyuchenkov & Shustov 
12004]) and the dust continuum radiative transfer code MC3D 
( jWolht al.|199^|Wolf|2003] ). 


Appendix B.3.1: URAN(IA) 


URAN(IA) is a 2D radiative transfer code for molecular lines. It 
features a Monte Carlo algorithm for calculating the mean inten¬ 
sity and an accelerated lambda iteration (ALI) method for self- 
consistent molecular excitation calculations. This code has been 
extensively tested in ID and 2D, for example against the RA- 
TRAN code ( jHogerheijde & van der Tak|2000| ). It successfully 
passed all tests formulated in t he benchmark paper for N-L TE 
molecular radiative transfer by jvan Zadelhoff et ak] p002| ). It 
has been used in several applications, for example modeling the 
starless core LI544 ( jPavlyuchenkov & Shustov|[2004| ), and in 
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- U RAN (I A) ART 

- U RAN (I A) LTE 

+ + Mol3D LTE 
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+ + Mol3D FEP 
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Fig. B.l. HCO^ (4-3) transition of a protoplanetary disk, oriented face- 
on. Different colors are obtained with different methods. Solid lines 
are obtained using the URAN(IA) software package and dotted lines 
represent solutions obtained with Mol3D. Shown is the normalized flux 
over velocity. Both codes produce comparable results within acceptable 
differences (less than 0.5%). LVG is a good appro ximation in the case 
of circumstellar disks ( [Pavlyuchenkov et al.|2007| ). 


dedicated parameter space studies of molecula r line formation 
of molecular line formation of prestellar cores (|Pavlyuchenkov| 
let al.|2QQ8D . 

In this code, several approximate methods of calculating the 
level populations are included. In a comprehensive study, 
jPavlyuchenkov et aL] ( |2QQ7| ) have shown the reliability of these 
methods in p rotoplanetary disk scenarios. 

In Sect. |2.1| we introduced the different methods of calculating 
level populations implemented in Mol3D. In Fig. |B.l|we present 


the spectrum of the HCO^ (4-3) transition of a typical protoplan¬ 
etary disk (Mdisk = 0-02 M©), comparing the different methods 
assuming a uniform N(HCO^)/N(H) ratio of 1-10“^. For sim¬ 
plicity, we choose a flared disk model around a T Tauri star uti¬ 
lizing a temperature distribution with a gradient in radial and 
vertical direction and a power-law model for the density distri¬ 
bution adapted from |Shakura & Sunyaev| ( p77^ : 


Pdustier) ~ ^ 

hir) = lOAU • ^ 


lOOAU 
r 


) expf 




h(ry 


(B.16) 


lOOAUJ 


The disk is orientated face-on. We assume pure Keplerian rota¬ 
tion and neglect dust re-emission. All spectra are normalized to 
the maximum value of the LTE solution. To allow comparison, 
the resulting spectra obtained with the URAN(IA) code are also 
included. 


URAN(IA) andMoUD produce comparable HCO^ (4-3) spectra 
with differences of less than about 0.5%. The most realistic so¬ 
lution (black line) is obtained using the accelerated Monte Carlo 
(ART), to classify the approximation methods. As discussed in 
greater detail in Appendix |B. 2. 1[ the LTE method overestimates 
the net flux and the EEP method underestimates it. The fluxes in 
the line obtained with the LVG method are in between those ob¬ 
tained with the other methods and are closely comparable to the 
ART method. This result has also been found by |Pavlyuchenkov 


et al.| ( |200'7] ) and we refer the interested reader to their study. 



Distance [AU] 


Fig. B.2. Midplane temperature of a typical protoplanetary disk 
around a pre-main-sequence T Tauri star. This disk inner radius 
amounts to 2 AU and the outer radius to 200 AU, respectively. The 
total disk mass amounts to 0.04 M©. Mol3D (red) and MC3D (blue) 
obtain comparable results. In this case, the differences depend signifi¬ 
cantly on the optical disk properties and subsequently on the number of 
photon packages simulated in the Monte Carlo process. Thus, the error 
is higher in the dense region at the inner rim of the disk (max. 10%) and 
lower at the outer disk regions (max. 2%). 


Appendix B.3.2: MC3D 

MC3D is a 3D continuum radiative transfer code that uses a 
Monte Carlo method to calculate self-consistent dust tempera¬ 
ture distributions, continuum re-emission/scattering maps, and 
SEDs. Among a wide variety of radiative transfer studies it 
has been applied extensively for the analysis of multiwavelength 
high-angular-resolution observations of circumstellar disks (e.g., 


Schegerer et al. 

20081 Sauter & WoltpOl ll|Madlener et al.|2012[ 

Grafe & Wolf] 

2013|). MC3D was successfully tested against 


other continuum RT codes (e.g., |Pascucci et al.|2004| ). 

In this section, we compare Mol3D with MC3D for a typical pro¬ 
toplanetary disk scenario. The density distribution is described 


by Eq. |B.16| assuming a total disk mass of 0.04 M©. The tem¬ 
perature structure is calculated self-consistently with each code 
using their Monte Carlo scheme. Both codes use a similar ap¬ 
proach based on the assumption of local thermal equilibrium and 
immediate temperature correction as proposed by |Bjorkman & 
|Wood| ( |200T] ). 

The resulting midplane temperature around a pre-main-sequence 
T Tauri star is shown in Eig. |B.2[ The calculation of the temper¬ 
ature in the dense midplane is one of the most crucial RT prob¬ 
lems, because these disk regions can hardly be reached by stellar 
photons directly. Thus, it is mostly heated indirectly by thermal 
re-emission radiation from the innermost disk regions and upper 
disk layers. Both codes produce very smooth and similar tem¬ 
perature distributions. The statistical nature of the Monte Carlo 
method results in deviations of about 10% for the innermost disk 
parts and less than 2% for the outer regions, which is mainly due 
to the probability that the photon packages can reach these re¬ 
gions. 

We also compared the dust re-emission/scattered light maps and 
SEDs and And that both codes produce comparable results with 
maximum deviations of about 10%. 
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